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LONG-TERM  GOALS 

The  long-term  goal  of  this  project  is  to  design  and  evaluate  the  components  that  will  comprise  a  next  genera¬ 
tion  mesoscale  atmospheric  model  within  the  Coupled  Ocean/Atmosphere  Mesoscale  Prediction  System 
(COAMPS®1 ').  It  is  anticipated  that  in  order  to  meet  future  Navy  requirements,  next  generation  approaches  to 
numerical  techniques  and  physical  parameterizations  will  be  needed. 

OBJECTIVES 

The  objectives  of  this  project  involve  the  development,  testing,  and  validation  of:  i)  new  numerical 
techniques  such  as  advection  schemes  and  time  differencing  methods,  and  ii)  new  methods  for  repre¬ 
senting  cloud-scale  physical  processes.  Both  of  these  objectives  are  tailored  to  address  high-resolution 
applications  for  horizontal  grid  increments  at  1  km  or  less. 

APPROACH 

Our  approach  is  to  follow  a  methodical  plan  in  the  development  and  testing  of  a  nonhydrostatic  micro¬ 
scale  modeling  system  that  will  leverage  the  existing  COAMPS  and  new  model  prototypes.  Our  work 
on  numerical  methods  will  involve  investigation  of  spatial  and  temporal  discretization  algorithms  that 
are  superior  to  the  current  generation  leap-frog,  second-order  accurate  numerical  techniques  presently 
employed  in  COAMPS  and  many  other  models;  these  new  discretization  methods  will  be  developed 
and  implemented.  Our  work  on  the  physics  for  the  next-generation  COAMPS  will  feature  the  devel¬ 
opment  of  physical  parameterizations  specifically  designed  to  represent  cloud-scale  processes  operat¬ 
ing  on  fine  scales.  A  parameterization  is  proposed  that  properly  represents  the  coupled  nature  between 
the  turbulence  and  microphysics  in  droplet  activation,  evaporation,  and  auto-conversion  processes  for 
mesoscale  and  microscale  models.  Validation  and  evaluation  of  the  modeling  system  will  be  per¬ 
formed  using  datasets  of  opportunity,  particularly  in  regions  of  Navy  significance. 

WORK  COMPLETED 

1.  Development  of  a  turbulence-based  cloud  droplet  activation  parameterization. 

We  developed  a  framework  for  parameterization  of  cloud  droplet  activation  in  a  high  resolution 
mesoscale  model.  There  are  two  key  elements  in  the  parameterization.  First,  the  activation  rate  is  a 
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strong  function  of  the  turbulence  intensity  represented  by  the  turbulence  probability  density  function 
(PDF).  Second,  the  activation  depends  on  the  activation  time  scale  that  is  defined  as  the  time  required 
for  a  parcel  to  reach  the  level  of  maximum  supersaturation  from  the  cloud  base.  This  time  scale  is  a 
function  of  vertical  velocity  and  CCN  spectrum.  Therefore,  the  new  activation  parameterization 
framework  is  closely  coupled  to  the  turbulence  structure.  We  included  this  parameterization  in  a  single 
column  high-resolution  turbulence  closure  model  and  tested  the  performance  against  DYCOMS  (Dy¬ 
namics  and  Chemistry  of  Marine  Stratocumulus  Cloud  Experiment  II)  observations. 

2.  Spectral  element  and  discontinuous  Galerkin  2D  prototypes. 

Spectral  element  (SE)  (see  Giraldo  2005)  and  discontinuous  Galerkin  (DG)  (see  Giraldo  2006)  meth¬ 
ods  are  a  new  class  of  spatial  discretization  methods  that  are  used  to  approximate  the  derivatives  of  the 
governing  equations.  In  COAMPS  and  WRF,  presently  this  is  done  using  the  finite  difference  method. 
The  advantage  of  SE  and  DG  methods  is  that  they  offer  high-order  accuracy  (unlike  low  order  finite 
differences)  and  this  accuracy  can  be  achieved  on  any  unstructured  grid  -  this  is  not  true  for  finite  dif¬ 
ference  methods  where  the  differencing  stencils  require  a  certain  level  of  structure  (such  as  orthogonal¬ 
ity).  Based  on  SE  and  DG  methods,  we  developed  5  different  prototypes  for  the  nonhydrostatic  Euler 
equations.  The  goal  of  this  exhaustive  study  was  to  detennine  not  only  which  method  to  use  for  a  next- 
generation  model  but  also  to  determine  which  set  of  equations  should  be  used.  For  example,  currently 
in  the  literature  at  least  three  different  fonns  of  the  Euler  equations  can  be  found.  The  first  set  we  de¬ 
note  as  set  number  1  and  is  written  as  follows: 

—  R  — 
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where  tenns  with  a  bar  denote  vector  quantities.  Equation  set  1  is  the  form  currently  used  in  many 
mesoscale  models,  including  COAMPS.  This  set  is  quite  popular  because  it  is  completely  self- 
contained;  meaning  that  there  are  three  equations  and  three  independent  variables  with  no  equation  of 
state  required  to  close  the  system.  The  solution  variables  in  this  case  are  exner  pressure  (tc),  velocity 
(u),  and  potential  temperature  (9).  This  equation  set  does  not  formally  conserve  mass  and  for  this  rea¬ 
son  it  makes  little  sense  to  either  write  the  equations  in  flux-form  or  to  use  conservative  methods  to 
solve  them.  For  example,  this  equation  set  is  ideally  suited  for  either  finite  difference  or  finite  element 
methods  which  conserve  mass  only  in  a  global  sense. 

A  set  that  is  becoming  increasingly  popular  is  set  number  2  which  is  written  as: 
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where  solution  variables  are  density  (p),  momentum  (pu),  and  potential  temperature  (p9);  however,  in 
this  form  of  the  Euler  equations  there  are  three  equations  and  4  unknown  variables  (p,  u,  0,  and  pres¬ 
sure).  Thus  in  this  case  an  equation  of  state  is  required  to  close  the  system.  It  is  important  to  note  that 
this  form  of  the  Euler  equations  is  written  in  flux-form  (conservation  form)  and  that  the  variables  are  in 
fact  conserved  quantities.  Therefore,  for  this  equation  set  it  makes  perfect  sense  to  use  locally  conser¬ 
vative  methods  such  as  finite  volume  or  discontinuous  Galerkin  methods. 

The  equation  set  number  3  studied  is  the  following: 

f-v.M 
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where  the  conservation  variables  are  density  (p),  momentum  (pu),  and  total  energy  (pe);  this  is  the  set 
that  is  typically  used  in  computational  fluid  dynamics  (CFD).  This  set  is  also  written  in  conservation 
form  where  all  the  variables  are  conserved  quantities.  While  this  set  has  been  very  popular  in  CFD,  it 
is  not  so  useful  for  geophysical  fluid  dynamics  since  most  parameterization  packages  are  written  in 
tenns  of  (virtual)  potential  temperature  and  not  total  energy.  While  we  developed  very  accurate  mod¬ 
els  using  this  equation  set,  we  eliminate  it  from  further  consideration  due  to  the  difficulty  of  having  to 
convert  from  total  energy  to  potential  temperature  between  the  dry  dynamics  and  moist  physics. 

3.  Weighted  Essentially  Non-Oscillatory  (WENO)  methods. 

Atmospheric  models  require  numerical  methods  that  can  accurately  represent  the  transport  of  tracers 
with  steep  gradients,  such  as  those  that  occur  at  cloud  boundaries  or  the  edges  of  chemical  plumes.  In 
atmospheric  sciences,  the  most  widely  used  numerical  techniques  for  this  type  of  problem  are  flux- 
corrected  transport  or  closely  related  flux-limiter  methods.  The  limiters  are  typically  designed  to  pre¬ 
vent  the  development  of  new  extrema  in  the  concentration  field.  This  will  preserve  the  non-negativity 
of  initially  non-negative  fields,  which  is  essential  for  the  correct  simulation  of  cloud  microphysics  or 
chemical  reactions.  One  serious  systematic  weakness  of  flux  limiter  methods  is  that  they  also  tend  to 
damp  the  amplitude  of  extrema  in  smooth  regions  of  the  flow,  such  as  the  trough  of  a  well-resolved 
sine  wave.  To  avoid  this  problem,  we  have  been  investigating  the  application  of  WENO  (Weighted 
Essentially  Non-Oscillatory)  methods  to  tracer  transport  in  atmospheric  models.  WENO  methods  are 
widely  used  in  many  disciplines,  but  scarcely  been  tested  in  atmospheric  applications.  WENO  meth¬ 
ods  preserve  steep  gradients  while  simultaneously  avoiding  the  dissipation  of  smooth  extrema  by  esti- 


mating  the  value  of  the  solution  in  a  way  that  heavily  weights  the  smoothest  possible  cubic  polynomial 
fit  to  the  local  function  values.  Where  the  solution  is  well  resolved,  all  possible  cubic  interpolants  are 
weighted  almost  equally.  Near  a  steep  gradient,  those  interpolants  that  straddle  the  gradient  are  almost 
completely  ignored. 

RESULTS 

1.  Development  of  a  turbulence-based  cloud  droplet  activation  parameterization. 

Figure  la  shows  the  dependence  of  the  activation  time  scale  on  the  vertical  velocity  and  CCN  spec¬ 
trum.  It  ranges  from  1.3  minutes  for  the  marine  spectrum  and  weak  updrafts  to  just  8  seconds  for  the 
continental  and  strong  updrafts.  The  time  scale  is  shorter  for  a  stronger  upward  motion  because  it  pro¬ 
vides  more  adiabatic  cooling;  and  it  is  longer  in  a  marine  air  environment  (lower  CCN  number)  since 
the  supersaturation  level  is  higher  and  the  condensation  rate  is  lower  than  those  in  continental  air  mass 
(higher  CCN  number).  The  turbulence  closure  model  with  the  activation  parameterization  is  used  to 
simulate  the  vertical  distribution  of  cloud  droplet  number  and  mean  radius  observed  in  DYCOMS2 
field  experiment.  The  modeled  droplet  number  concentration  and  turbulence  variables  agree  well  with 
the  observations.  Our  simulations  suggest  that  the  new  parameterization  should  be  applicable  to  next 
generation  high  resolution  COAMPS  where  turbulence  PDF  may  be  well  predicted  or  diagnosed. 


Activation  Time  Scale 


W  (cm/s) 


Droplet  Number  Density 


Figure  1.  (a):  Activation  time  scale  as  a  function  of  upward  velocity  and  CCN  spectrum.  The 
different  CCN  spectrum  used  in  the  calculation  is  denoted  by  “ Marine ”,  " Continental ”,  and  “C-M” 
(the  mixed  air  mass)  respectively;  (b):  Vertical  profile  of  simulated  cloud  droplet  number  density. 
The  circles  are  leg  averages  of  aircraft  data  taken  in  DYCOMS  flight  1;  and  horizontal  bar  is  the 

scatter  of  the  data. 


2.  Spectral  element  and  discontinuous  Galerkin  2D  prototypes. 


In  order  to  determine  which  equation  set  along  with  which  numerical  method  is  best  suited  for  building 
a  next-generation  mesoscale  model,  we  ran  simulations  for  two  standard  test  cases.  The  test  cases  se¬ 
lected  were  the  rising  thermal  bubble  problem  made  famous  by  Andre  Robert  and  the  linear  hydro¬ 
static  mountain  wave  used  by  Durran,  Klemp,  Skamarock,  and  analyzed  in  detail  by  Ron  Smith.  We 
used  the  spectral  element  and  discontinuous  Galerkin  methods  to  solve  the  equations.  We  only  show 
results  for  the  SE  method  because  both  methods  yield  very  similar  results  for  these  two  particular  test 


cases.  However,  it  is  expected  that  for  more  challenging  problems  (i.e.,  those  having  very  steep  gradi¬ 
ents)  the  DG  method  should  prevail. 

The  time-integrators  used  to  advanced  the  models  are  a  family  of  strongly  stability  preserving  (SSP) 
schemes  which  are  sometimes  referred  to  as  TVD  time-integrators  (TVD=  total  variation  diminishing). 
These  time-integrators  control  unwanted  spurious  oscillations  near  large  gradients  and  when  used  in 
combination  with  slope  limiters  with  the  DG  method  results  in  a  truly  TVD  method  in  both  space  and 
time.  The  TVD  property  is  important  because  it  means  that  the  models  never  produce  unphysical  ex¬ 
trema  throughout  the  time-integration  regardless  of  the  strength  of  the  gradients.  These  methods  have 
not  been  used  previously  in  atmospheric  models  and  we  use  a  family  of  2nd  and  3ld  accurate  SSP  meth¬ 
ods. 

2a.  Rising  thermal  bubble. 

The  rising  thermal  bubble  problem  does  not  have  an  analytic  solution  since  it  involves  the  full  nonlin¬ 
ear  equations  but  is  nonetheless  a  useful  test  because  the  solution  is  intuitively  straightforward  to  un¬ 
derstand.  Another  attractive  trait  of  this  problem  is  that  the  boundary  conditions  are  quite  simple  since 
they  only  require  no-flux  across  the  domain  boundaries  which  gives  a  good  measure  of  how  the  dis¬ 
crete  operators  are  behaving. 

Figure  2  shows  the  color  contours  for  the  potential  temperature  perturbation  from  the  isothermal  refer¬ 
ence  state  after  600  seconds.  The  result  shown  was  obtained  with  the  spectral  element  model  using 
equation  set  one;  however,  this  result  is  identical  for  all  three  equation  sets  using  either  the  SE  or  the 
DG  method.  The  initial  thermal  perturbation  is  a  cosine  wave  with  maximum  of  0.5  above  the  refer¬ 
ence  state.  Note  that  the  color  contours  seen  in  Figure  2  are  between  0  to  0.5  which  shows  that  the 
model  does  not  produce  spurious  extrema  (overshoots  or  undershoots).  In  addition,  the  solution  is  very 
symmetric  about  the  x=500  meter  axis  which  is  what  one  would  expect. 


a)  solution  at  300  seconds  b)  solution  at  600  seconds  c)  solution  at  900  seconds 


Figure  2:  The  potential  temperature  perturbation  from  the  isothermal  reference  state  after  a)  300, 
b)  600,  and  c)  900  seconds  for  the  spectral  element  and  discontinuous  Galerkin  nonhydrostatic 
models  using  a  grid  of  160x160  grid  points  in  a  1  km  x  1  km  domain. 


2b.  Linear  hydrostatic  mountain  waves. 

The  linear  hydrostatic  mountain  wave  problem,  on  the  other  hand,  has  an  analytic  solution  which  al¬ 
lows  us  to  discuss  quantitatively  the  performance  of  specific  models.  The  main  issue  with  this  test  case 


is  that  it  requires  more  sophisticated  boundary  conditions  (such  as  either  non-reflecting/radiative 
boundary  conditions  or  sponge  layers)  in  order  to  properly  perform  the  simulation.  In  the  end,  the  ac¬ 
curacy  of  the  simulation  is  completely  determined  by  the  quality  of  the  lateral  and  top  boundary  condi¬ 
tions.  Sophisticated  radiation  boundary  conditions,  while  useful  for  this  test  case,  unfortunately  are 
not  so  useful  in  an  operational  setting  where  a  global  model  is  used  to  drive  the  boundary  conditions  of 
the  mesoscale  model. 


b)  Analytic  Solution 


Figure  3:  The  vertical  velocity  for  the  linear  hydrostatic  mountain  wave  after  10  hours  for  a)  the 
spectral  element  model  and  b)  the  linear  hydrostatic  analytic  solution.  The  simulation  uses  160  x 

160  grid  points  for  a  240  km  x  30  km  domain. 

Figure  3  shows  the  contours  for  the  vertical  velocity  after  a  10  hour  simulation  for  the  linear  hydro¬ 
static  mountain  with  a  height  of  one  meter  and  half  width  of  10  kilometers.  The  numerical  solution  is 
shown  on  the  left  panel  and  the  analytic  solution  on  the  right. 


Defining  the  root-mean-square  (RMS)  error  as  follows: 


RMS 


where  the  superscripts  N  and  A  denote  the  numerical  and  analytic  solutions,  and  Np  is  the  number  of 
grid  points,  we  compute  the  RMS  error  for  the  vertical  velocity  to  be  1.09xl0‘4  which  is  extremely 
competitive  with  the  results  obtained  with  WRF  and  COAMPS. 


3.  Weighted  Essentially  Non-Oscillatory  (WENO)  methods.. 


The  following  three  figures  illustrate  how  WENO  methods  can  perform  in  an  atmospheric  context  in¬ 
volving  trapped  mountain  lee  waves.  Fig.  4  shows  the  isentropes  of  potential  temperature  and  the  ver¬ 
tical  velocity  field  (color  fill)  in  simulations  using  numerical  models  that  differ  only  in  their  treatment 
of  potential  temperature  and  passive  tracer  transport.  In  the  top  panel  (a)  ,  the  WENO  method  is  used 
for  the  advection,  in  the  middle  panel  (b)  leapfrog  time,  centered  4th-order  spatial  difference  is  used, 
and  in  the  bottom  panel  (c)  a  flux  limiter  method  proposed  by  LeVeque  is  employed.  The  leapfrog 
scheme  is  non-damping  and  correctly  produces  a  large  amplitude  lee-wave.  The  WENO  method  gives 


essentially  the  same  solution  as  the  leapfrog  scheme,  but  the  flux-limiter  method  is  incorrectly 
damped. 

Now  consider  a  complimentary  test  in  which  a  passive  tracer  with  uniform  concentration  of  1 .0  is  ini¬ 
tially  distributed  throughout  a  rectangular  region  upstream  of  the  waves.  There  are  sharp  gradients  at 
each  edge  of  the  distribution,  as  shown  in  Fig.  5.  As  this  tracer  advects  downstream,  it  is  distorted  by 
the  shear  and  curvature  in  the  wind  field,  but  its  edges  should  remain  sharp  and  it  should  remain  be¬ 
tween  the  two  isentropes  that  initially  passed  above  and  below  it.  As  apparent  in  Fig.  6,  the  leapfrog 
scheme  (panel  b)  is  not  capable  of  handling  the  steep  gradients  in  the  tracer  distribution  and  creates  a 
series  of  undershoots  and  overshoots.  On  the  other  hand  both  the  WENO  and  the  flux-limited  methods 
perform  well,  with  the  WENO  method  preserving  steeper  gradients  than  the  flux  limiter  method. 

These  tests  show  that  the  WENO  method  is  capable  of  performing  well  in  both  situations  (smooth  flow 
with  important  extrema,  and  tracer  transport  with  steep  gradients),  whereas  each  of  the  other  methods 
only  works  well  one  of  the  cases. 


Figure  4.  Isentropes  and  the  vertical  velocity  field  (color  fill)  for  (a)  WENO 
method  is  used  for  the  advection,  (b)  leapfrog  time,  centered  4th-order  spatial  dif¬ 
ference,  and  (c)  a  flux  limiter  method  proposed  by  LeVeque. 
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Figure  5.  Initial  passive  tracer  distribution  with  uniform  concentration  of  1.0  as 
represented  by  the  rectangular  region  upstream  of  the  waves. 


Figure  6.  Isentropes  and  the  passive  tracer  (color  fill)  for  (a)  WENO  method  is 
used  for  the  advection,  (b)  leapfrog  time,  centered  4th-order  spatial  difference,  and 
(c)  a  flux  limiter  method  proposed  by  LeVeque. 


IMPACT/APPLICATIONS 


COAMPS  is  the  Navy’s  operational  mesoscale  NWP  system  and  is  recognized  as  the  key  model  com¬ 
ponent  driving  a  variety  of  DoD  tactical  decision  aids.  Accurate  mesoscale  prediction  is  considered  an 
indispensable  capability  for  defense  and  civilian  applications.  Skillful  COAMPS  predictions  at  resolu¬ 
tions  less  than  1  km  will  establish  new  capabilities  for  the  support  of  the  warfighter  and  Sea  Power  2 1 . 
Operational  difficulties  with  weapon  systems  such  as  the  Joint  Standoff  Weapon  (JSOW)  have  been 
documented  in  regions  with  fine-scale  topography  due  to  low-level  wind  shear  and  turbulence.  Im¬ 
proved  high-resolution  predictive  capabilities  will  help  to  mitigate  these  problems  and  introduce  poten¬ 
tially  significant  cost  saving  measures  for  the  operational  application  of  JSOW.  The  capability  to  pre¬ 
dict  the  atmosphere  at  very  high  resolution  will  further  the  Navy  sea  strike  and  sea  shield  operations, 
provide  improved  representation  of  aerosol  transport,  and  will  lead  to  tactical  model  improvements. 
Applications  of  CO  AMPS  at  resolutions  less  than  1  km  will  establish  important  direction  for  the  de¬ 
velopment  of  the  Navy's  next  generation  microscale  prediction  system.  Emergency  response  capabili¬ 
ties  and  Homeland  Security  issues  within  the  DoD  and  elsewhere,  such  as  LLNL,  will  be  enhanced 
with  the  new  modeling  capability. 

TRANSITIONS 

The  next  generation  COAMPS  system  will  transition  to  6.4  projects  within  PE  0603207N  (SPA WAR, 
PMW-180)  that  focus  on  the  transition  COAMPS  to  FNMOC. 

RELATED  PROJECTS 

COAMPS  will  be  used  in  related  6. 1  projects  within  PE  0601 153N  that  include  studies  of  air-ocean 
coupling,  boundary  layer  studies,  and  topographic  flows  and  in  related  6.2  projects  within  PE 
0602435N  that  focus  on  the  development  of  the  atmospheric  components  (QC,  analysis,  initialization, 
and  forecast  model)  of  COAMPS.  . 
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